Spin-Peierls transition in the Heisenberg chain with finite-frequency phonons 
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We study the spin-Peierls transition in a Heisenberg spin chain coupled to optical bond-phonons. 
Quantum Monte Carlo results for systems with up to N — 256 spins show unambiguously that the 
transition occurs only when the spin-phonon coupling a exceeds a critical value a c . Using sum rules, 
we show that the phonon spectral function has divergent (for N — > oo) weight extending to zero 
frequency for a < a c . The equal-time phonon-phonon correlations decay with distance r as 1/r. 
This behavior is characteristic for all < a < a c and the q — it phonon does not soften (to zero 
frequency) at the transition. 

PACS numbers: 75.10.Jm, 75.40. Cx, 75.40.Mg, 63.22.+m 
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The S — 1/2 Heisenberg spin chain is unstable towards 
dimcrization (the spin-Peierls transition) when coupled 
to an clastic lattice H. For phonons in the adiabatic 
limit, this transition has been predicted to occur for arbi- 
trarily weak spin-lattice coupling |L| . On the other hand, 
recent work for optical phonons in the anti-adiabatic 
(high-frequency) limit suggests a transition only above 
a critical coupling J2|. Furthermore, the mechanism of 
the transition in this limit was suggested to be qualita- 
tively different, with no softening of the q — n phonon 
g. A way to reconcile the results in the adiabatic and 
anti-adiabatic limits has been proposed within an im- 
proved mean- field (RPA) theory ||, with the result that 
a complete phonon softening occurs only for bare phonon 
frequencies loq less than a critical value. For higher o>o, 
a central peak appears in the phonon spectral function 
and the phonon branch remains gapped. Considering the 
manifestly uncontrolled nature of mean-field calculations 
in one dimension, non-perturbative results in the regime 
of phonon frequencies comparable to the magnetic ex- 
change energy J are required to test this novel scenario. 

In this Letter, we address the issues of a critical 
spin-phonon coupling and the mechanism of the zero 
temperature dimerization transition in the strictly one- 
dimensional case, using quantum Monte Carlo (QMC) 
simulations to obtain numerically exact results for rela- 
tively large systems. The model we study is defined by 
the Hamiltonian 
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is the phonon coordinate, and rij = ajcii is the phonon 
occupation number at bond i. We use a recently devel- 
oped QMC method based on sampling the perturbation 
expansion in the interaction representation B. For a fi- 
nite lattice at finite inverse temperature /3, the expansion 



converges for any decomposition of H — Hq + V into di- 
agonal (Ho) and perturbing (V) terms and can be used 
H as a basis for a "worldline" Monte Carlo algorithm 
in continuous imaginary time (i.e., without invoking the 
Trotter decomposition |J) . In a slight modification of the 
scheme introduced in Ref. g, we here include only the 
bare phonons in the diagonal term; H — uj J2i n i- The 
updating of the spin degrees of freedom can then be car- 
ried out using a new and highly efficient "operator-loop" 
algorithm |7j], which in particular allows for sampling of 
all winding number sectors and hence direct evaluation 
of the spin stiffness || . 

In this work we consider only an energy ujq = J/4 for 
the bare phonons and study the behavior for values of the 
spin-phonon coupling in the range < a/ J < 0.5. We 
have studied systems with N up to 256 at inverse temper- 
atures (3 = J/T sufficiently high to give ground state re- 
sults. Typically, for the system sizes we have considered, 
P as high as rs 2A^ is required to achieve convergence to 
the T = limit of all the quantities of interest. We have 
used at least j3 = AN for all calculations presented here. 

Note that in the model, Eq. p), for a > there 
is an energy gain associated with an average uniform 
phonon displacement (x) = (l/N)J2i( x i) > 0, which 
leads to an increased average effective spin-spin coupling 
J cS = J + a (x) > J. For cj /J = 0.25 at T = 0, 
we find J off / J = 1.018,1.071,1.158,1.278, and 1.430 for 
a = 0.1,0.2,0.3,0.4, and 0.5, respectively. We will in 
some cases measure energies in units of J c ff instead of 
the bare exchange J. 

The most direct signal of the dimerization that can be 
measured in our simulations is the approach of the stag- 
gered phonon-phonon correlation function (— l) r (xiXi+ r ) 
to a non-zero value at long distances r. Our results for 
this quantity indicate a critical coupling 0.1 < a c /J < 
0.35 for ujq/J = 0.25 ||. In order to improve on the 
accuracy of this rough estimate, and to circumvent po- 
tential problems with detecting a very small dimcrization 
(as would be the case for very weak coupling if the mean- 
field result a c — would be correct), we have also con- 



0.30 



0.25 - 



0.20 



0.15 - 



0.10 



0.05 



0.00 
0.00 




0.10 



0.20 



0.30 



0.40 



0.50 



a/J 



FIG. 1. Spin stiffness vs the spin-phonon coupling for sys- 
tem sizes N — 8, 16, 32, 64 and 128 (p s decreases with increas- 
ing N). Where not shown, statistical errors are smaller than 
the symbols. The curves are high-order polynomial fits to the 
data points. 

sidered several other quantities. It is particularly useful 
to study the effects of the dynamic phonons in the spin 
sector (in general, our simulation results for spin quanti- 
ties have smaller statistical fluctuations than the phonon 
correlations). Thus we discuss here results for the spin 
stiffness and the staggered spin susceptibility. 

If the dimerization transition occurs at some critical 
coupling a c > 0, it is expected to be of the Kosterliz- 
Thouless (KT) type JlCfl , The spin stiffness p s , i.e., the 
ground state energy curvature with respect to a uniform 
twist <j) in the spin-spin interaction pT| , 
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is then expected to exhibit a discontinuous jump from 
a finite value for a < a c to zero for a > a c (reflect- 
ing the opening of a spin gap). For a finite system the 
jump will be smoothed. In Figure \u we show results for 
the stiffness versus a/J for several system sizes. The 
behavior expected for a KT transition is seen clearly — 
p s rapidly approaches zero for a/J > 0.4 but appears 
to converge to a finite value for a/J < 0.2, indicating a 
critical coupling between these values, in agreement with 
our previous results for the dimerization. It is, however, 
not easy to extract an accurate value for a c using these 
results. The scaling behavior is complicated by logarith- 
mic corrections, which we expect to be present for all 
a < a c as in the case of the Heisenberg chain [i.e., a = 
in Eq. O)]. This is in contrast to the finite tempera- 
ture KT transition in the two-dimensional XY model, 
where p s approaches its asymptotic value algebraically 
for T < T c and logarithmically only exactly at T c |13,E3[ ■ 
An indication of the difficulties associated with the log 
corrections in the spin-phonon chain can be seen in our 
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FIG. 2. Size dependence of the staggered spin susceptibil- 
ity for different values of the spin-phonon coupling. A linear 
behavior of [x s {^)/N] 2 vs \n(N) is expected in the gap less 
phase, with slope zero at the critical point. The decrease 
with increasing N for a/J = 0.3 indicates the presence of a 
spin gap. 



stiffness data for a = 0, for which the exact infinite-size 
value is known to be (in our units) pTI] p s — 1/4; about 
8% lower than what we find for N = 128. 

Although the log corrections complicate the extraction 
of a c from the stiffness data, their presence in other quan- 
tities can in fact be useful in numerical calculations. The 
asymptotic behavior of the spin-spin correlation function 
of the Heisenberg chain is known from bosonization and 
conformal field theory |lj.|l5[ ; 
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m 1/2 (r/r ). 



(4) 



We expect this form to apply for all a < a c . The loga- 
rithmic correction should vanish at the critical point a c , 
as it is known to do, e.g., at the critical point of the 
frustrated Ji — J 2 chain 16 
staggered spin susceptibility 



We have calculated the 
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for which Eq. (|4|) and conformal invariance imply the 
finite-size scaling form |1J] 



X a (7r)~-Mn 1/2 (JV/iVo). 
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In Figure @ we graph (x s (-7r)/./V) 2 vs ln(A) for a/J in 
the range 0.1 — 0.3. For a/J = 0.1 and 0.2 the linear 
behavior for the larger system sizes is consistent with the 
form (0) expected in the gapless phase, whereas in the 
a/J = 0.3 case there is a clear decrease with increasing 
A, corresponding to a finite asymptotic value for Xs{^) 
and therefore the presence of a spin gap. For a/J — 0.23 



the curve is flat within statistical errors for N > 64, 
implying that Xs(^) diverges linearly with N without log 
correction. Based on this result (and calculations for 
other a/ J close to 0.23) we conclude that a c /J = 0.23 ± 
0.01. 

Having established a KT transition and the critical 
coupling, we now turn to the question of the behavior of 
the q = 7r phonons at the transition. We consider the 
phonon spectral function 

A(q,u) = J2^ En \(™\x q \n)\ 2 5(u - [E m - E n \), (7) 



where 
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This real-frequency dynamic quantity cannot be obtained 
directly in our simulations. In order to avoid the prob- 
lems associated with numerically continuing imaginary 
time data to real frequency, we here study sum rules that 
relate A(q, u>) to quantities that can be directly calcu- 
lated. Two useful integrals that can be easily obtained 
from Eq. ([fl) are 
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where S x (q) and Xx (q) are the static structure factor and 
susceptibility; 



S x (q)= (x^ q X q ) 

Xx(q)= / dT(x_ q (T)x q (0)). 
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Using Eqs. (pal) and (ptl) one can readily verify |17j that 
the ratio 



R(q) = 2S x (q)/ Xx (q) 
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is an upper bound for the lowest phononic excitation of 
momentum q. For a > a c we therefore expect R(ir) — + 
as N — > oo, reflecting the presence of two degenerate 
ground states with momenta and 7r (linear combina- 
tions of the two possible real-space dimerized states). For 
a = 0, R(q) = loq for all q. A transition caused by a soft- 
ening of the q — 7r phonon would imply R(n) > for 
a < a c and R(jr) — > as a — > a c . This behavior is 
not seen in our results. Instead, R(tt) appears to ap- 
proach zero as N is increased even for a much smaller 
than the critical value, as shown in Figure k5| Hence the 
spectral weight extends to zero frequency also in the non- 
dimerized systems. 

We also find that the total spectral weight, given by 
the structure factor according to Eq. (Bal), diverges with 
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FIG. 3. Upper bound R(n) — 2S x (n)/xx(Tr) for the lowest 

q = 7r phonon excitation energy vs the inverse system size. 

The dashed line is at the non-interacting (a = 0) value ljq/J. 



N. In Figure we graph S x (ir) versus the logarithm of 
the system size for values of a both below and above the 
critical coupling. We find a linear increase for a < a c , in- 
dicating a logarithmic divergence and therefore an inverse 
distance decay of the real-space phonon-phonon correla- 
tion function. For a = 0.3 > a c the divergence is faster 
than logarithmic. The expected linear in N behavior 
cannot be observed close to ct c for the system sizes we 
have studied, due to large short-distance contributions 
to S x (tt). For a > 0.4 we do observe an almost linear 
divergence with N. 

These results show that, in the thermodynamic limit, 
there is infinite q = -k phonon spectral weight also for 
a < a c . This weight extends to zero frequency. The rate 
of decay of R(tt) with increasing N, seen in Fig. y, shows 
that the low-frequency weight grows rapidly with N. The 
only plausible explanation for this is that the phonon 
spectral function has a central peak with infinite integral. 
We now elaborate on the reasons for this behavior. 

In the absence of spin-phonon couplings, the low- lying 
excitations of the system are the two-spinon singlet and 
triplet states of the Heisenberg chain. Our results indi- 
cate that an arbitrarily weak coupling to the phonons in- 
duces an infinite phonon spectral weight into these states 
at the staggered momentum q — ir. The asymptotic real- 
space staggered phonon correlation function has the same 
1/r decay as the spin-spin correlation function (perhaps 
differing by multiplicative logarithmic corrections that 
cannot be detected in our results for the phonon corre- 
lations). This is not completely surprising, considering 
that the Heisenberg chain is also characterized by an in- 
verse distance decay of the dimerization correlation func- 
tion ((Si ■ Si+i)(Si+ r ■ Sj+i+r)}, as also recently noted by 
Gros and Werner p[. The corresponding susceptibility 
is therefore divergent and this leads to the spontaneous 
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FIG. 4. Staggered phonon structure factor vs the logarithm 
of the system size. The linear behavior for a = 0.1 and 0.2 cor- 
responds to an inverse distance decay of the real-space phonon 
correlation function. Note that with our definition [Eq. (J2J)] 
of the phonon coordinate, the non- interacting structure factor 
equals 0.5. 



trix renormalization group study of a model closely re- 
lated to the one we have studied also finds unambigu- 
ously that the critical coupling a c > for any ujq p9[ . 
In addition, our results are consistent with the calcula- 
tions in the anti-adiabatic limit H (although the infinite 
low-frequency q — n phonon weight for a < a c was not 
noted there), even though our ojq = J/4 is closer to the 
adiabatic regime. 

We have here discussed only the T — quantum phase 
transition in the strictly one-dimensional case. The fi- 
nite T c in real materials such as CuGeC>3 [Q can be due 
to three-dimensional phonons, as well as interchain mag- 
netic couplings. The non-softening nature of the quan- 
tum phase transition that we have found here clearly 
supports the suggestion |g|| that the finite-T transition 
may also be non-softening. In the improved RPA the- 
ory [|3| , softening occurs below a critical value of ujq . In 
future work, we plan to extend our simulations to two- 
and three-dimensional systems and hope to address this 
important issue. 

We would like to thank Robert Bursill, Valeri Kotov, 
Ross McKenzie, Oleg Sushkov, and Johannes Voit for 
useful discussions. The numerical simulations were car- 
ried out at the NCSA facilities at the University of Illisois 
at Urbana- Champaign. This work is supported by the 
National Science foundation under Grant No. DMR-97- 
12765. 



dimerization for arbitrarily weak spin-phonon couplings 
in the adiabatic case ojq = 0. What we have shown here 
is that dynamic (loo > 0) phonons destroy the long-range 
order for weak spin-phonon couplings but nevertheless 
the spin and phonon excitations remain coupled in a man- 
ifestly non-perturbative fashion. 

The most plausible scenario for the mechanism of the 
spin-Peierls transition is then the following: For weak 
coupling a, the phonon spectral function at q = tt has 
a finite- weight peak close to the bare frequency ujq, as 
well as a central peak with divergent frequency integral 
[In (iV) divergent as a function of system size] . As a is 
increased the finite-frequency peak may shift slightly but 
remains at finite frequency. The central peak sharpens 
and at a = a c acquires a S function component, corre- 
sponding to the development of static long-range order. 
For a > a c this 6 function remains, and a gap develops 
to the remainder of what was the finite-width central 
peak for a < a c . This gap is the excitation energy of a 
lattice/magnetic soliton pair. 

On general grounds, we find it unlikely that there 
would be any qualitative changes in the nature of the 
T = transition as ujq is varied. A recent density ma- 
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